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ABSTRACT 

We perform spatial and spectral analyses of the LMC gamma-ray emission collected over 66 months 
by the Fermi Gamma-ray Space Telescope. In our spatial analysis, we model the LMC cosmic-ray 
distribution and gamma-ray production using observed maps of the LMC interstellar medium, star 
formation history, interstellar radiation field, and synchrotron emission. We use bootstrapping of 
the data to quantify the robustness of spatial model performance. We model the LMC gamma-ray 
spectrum using fitting functions derived from the physics of tt^ decay, bremsstrahlung, and inverse 
Compton scattering. We find the integrated gamma-ray fiux of the LMC from 200 MeV to 20 GeV to 
be 1.37±0.02 X 10 ^ ph cm ^ s of which we attribute about 6% to inverse Compton scattering and 
44% to bremsstrahlung. From our work, we conclude that the spectral index of the LMC cosmic-ray 
proton population is 2.4±0.2, and we find that cosmic-ray energy loss through gamma-ray production 
is concentrated within a few 100 pc of acceleration sites. Assuming cosmic-ray energy equipartition 
with magnetic fields, we estimate LMC cosmic rays encounter an average magnetic field strength ^3 
fiG. 


1. INTRODUCTION 

The LMC is by far the brightest source of diffuse 
gamma-ray flux detected outside of the Milky Way. The 
LMC was first observed in gamma rays by the Ener¬ 
getic Gamma Ray Experiment Tel escope with signifi¬ 
cance > 4.5cr ( Sreeknmar et al.|1992 ). However, a spatial 
comparison of the LMC gamma-ray emission with ob¬ 
servations at other wavelengths was not possible before 
the Fermi Gamma-ray Space Teles cope. The fi r st such 
spatial analysis was performed by Abdo et al. (2010a) 
who determined Fermi had detected the LMC at 33(7 
significance after 11 months of observation. As an exter¬ 
nal, star-forming galaxy detected at high significance in 
gamma rays, the LMC provides a cosmic-ray laboratory 
that complements the Milky Way, which has been exten¬ 
sively modeled by the Galactic cosmic-ray pro pagation 
team (GALPROP, |Strong & Moskalenko||1998|). 

The LMG is at a wel l -constrained distance of 50 kpc 
(jMatsunaga et al.||2009t iPietrzyhski et al.||2013|) with a 
modest inclination of approximat ely 30"" (|de Vaucouleurs 


fc Ereeman|1972[ Kim et al.|1998[[van der Marel ^ Cioni 


2001). This small inclination angle means self extinc- 
tion and source confusion are low in the LMG, and the 
galaxy’s high Galactic latitude means the LMG does 
not suffer greatly from extinction by the Milky Way. 
If unaccounted for, the LMG’s inclination angle intro¬ 
duces a <10% error on distance measurements, much less 
than distance uncertainties to some objects in the Milky 
Way, e.g., SNR RGW 86 for which various stu dies have 
adopted a distance anywhere from 2 to 3 kpc (|Morlino| 
et al.||2014 ). Eurthermore, the LMG is near enough that 
stars can be resolved and the interstellar medium (ISM) 
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can be studied with high resolution (a few pc), yet it is 
distant enough to allow a global view of its stars and ISM. 
The detailed and global knowledge of stars and ISM em¬ 
pirically anchor the analysis of gamma-ray observations 
and study of cosmic rays in the LMG. 

The LMG hosts active star formation, creating massive 
stars that live and explode as supernovae in a few to a 
few tens of Myr. The massive stars’ fast winds and final 
supernova explosions interact with the ambient ISM and 
generate shocks, w hich produce cosmic rays by diffusive 
shock acceleration (Eermi||1949 Ackermann et al.||2011 


Pavlidou & Eields (2001) predicted that the LMG 
should be the brightest external normal star-forming 
galaxy seen by Fermi, along with the SMG. Fermi has 
confirmed these predictions and firmly established star¬ 
forming galaxies as a gamma-ray source cla ss, which 
has lead to rec ent t heoretical stud ies such as IPersic fc 


Rephaelil (|2012|) and|Martin| (|2014|). The SMG (|Abdo et 
aTlSOTUbln^ (I Abdo et al. MUB, M82 and INGC 253 


([AbT 

[ Abdo et al.|2010d|) , and two other starburst galaxies are 
etected, but all of these sources do not have the spatial 
extent of the LMG and are detected at a much lower 
signal-to-noise ratio in gamma rays. 

The LMG stands out as the most intense, spatially ex¬ 
tended source of gamma-ray emission outside the Milky 
Way; thus, the Fermi gamma-ray observations of the 
LMG render it an invaluable laboratory for the study 
of cosmic-ray physics. Moreover, star-forming galaxies 
are important co ntributors to the extragalactic gamma- 


ray b ackground (Pavlidou & Eields 2002 Eields et al 


2010), and so an understanding of the gamma-ray emis¬ 


sion properties in individual resolved galaxies is impor¬ 
tant for modeling the unresolved contribution to the 
background. 

Gosmic rays propagate through their host galaxy until 
they either escape or lose their energy through photon 
emission and particle-particle collisions. Diffuse gamma- 
ray emission is often produced along the way by cosmic- 
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ray protons and electrons via three different mechanisms. 
Cosmic-ray protons lose energy via hadronic gamma-ray 
emission when they collide with thermal ISM protons to 
produce neutral pious that quickly decay to gamma rays: 

^’CR + PisM + 7r°, 7r° ^ 7 + 7. ( 1 ) 

Cosmic-ray electrons emit bremsstrahlung radiation at 
gamma-ray energies when they scatter off of interstellar 
nuclei: 

^CR + PiSM ^ e -1- -1- 7 (2) 

Cosmic-ray electrons can also up-scatter interstellar radi¬ 
ation photons to gamma-ray energies via inverse Comp¬ 
ton scattering: 


-CR 


+ 7ISRF ^ e -h 7 ^ — r 


(3) 


The spatial analysis performed by Abdo et ah (2010a) 
indicates the Fermi detection of the LMC is likely dii- 
fuse, as expected if these interactions are indeed the 
sources of the gamma-ray emission. The highly signif¬ 
icant, spatially extended Fermi detection of the LMC, 
complemented by detailed and global knowledge of the 
distribution of stars and ISM, allows us to empirically 
test the relative importance of these three mechanisms 
in generating gamma rays. 

Here we investigate the relative contributions to the 
LMC gamma-ray signal W each of the cosmic-ray in¬ 
teractions given in Eqs. 0-(§. We perform a binned 
likelihood analysis to determine where in the LMC these 
interactions occur. We expand on the spatial analysis of 
Abdo et al. (2010a) by: 


1. Using the distribution of LMC star formation and 
tracers of the LMC’s interstellar radiation field 
(ISRF) as quantitative probes of gamma-ray inten¬ 
sity. 


2. Enforcing stricter photon selection criteria for bet¬ 
ter angular resolution given the luxury of more 
data. 


3. Bootstrapping the Fermi data to test robustness of 
spatial model rankings. 

4. Applying, a cross the entire galaxy , the smoothing 
kernels used Murphy et al. (2012) in their cosmic- 
ray diffusion analysis of the LMC star-forming re¬ 
gion 30 Doradus (30 Dor). 

In an independent analysis, we assess the spatially 
integrated LMC gamma-ray spectrum as a means for 
predicting the flux contributions of the gamma-ray pro¬ 
ducing cosmic-ray interactions. We simultaneously fit 
three theoretical spectra (pion decay, bremsstrahlung, 
and inverse Compton) to the LMC spectrum. The fit- 
te d spectra are based on the on e-zone model presented 
by jChakraborty & Fieldsj (|2013|) in their study of inverse 
Compton gamma-ray emission in star-forming galaxies 
where we have rescaled those authors’ model of the Milky 
Way to appropriately model the LMC. Ou r method ex¬ 
pands on the work of Abdo et al. (2010a) by indepen¬ 
dently fitting each of tnese three spectra rather than as¬ 
suming their relative ratios and fitting a single spectral 
model of their sum. 


This article proceeds as follows: describes our data 

selection and preparation. In ^ we list our sources of 
background contamination and the means by which we 
remove them. Q addresses the possibility that a point 
source unassociated with the diffuse emission of the LMC 
may be in or behind the galaxy. In ^ we motivate 
our models for the LMC gamma rays. In ^ and ^ we 
present the results of our spatial and spectral analyses 
respectively. We discuss our results in the context of past 
studies and motivate future work in ^ and we make our 
concluding remarks in ^ 


2. DATA 


Following Abdo et al. ( |2QlQa ), we select a 20° x 20° re¬ 
gion surrounding the LMC (centered at 5^20™, —68°30'). 
This large region of interest that extends well beyond the 
LMC gamma-ray emission is crucial for modeling of back¬ 
ground sources (see ® . These data were collected by the 
Fermi Large Area T^scope from 2008 August 8 to 2014 
January 30, and the gamma rays we analyze range in en¬ 
ergy from 200 MeV to 20 GeV. Figure shows the raw 
data with a TAN projection (Calabretta & Greisen|2002 ). 
We chose this projection oi the celest i al sphe re, which 


differs from that used in Abdo et al. (2010a), for eas 


ier comparisons with images at other wavelengths, e.g. 
Figures and 



RA (J2000) 


Fig. 1. — Fermi data from a 20° X 20° region surrounding the 
LMC. Photons range in energy from 200 MeV to 20 GeV and were 
collected over a period of about 66 months. The dashed box shows 
the extent of the images in Figure!^ A square root color scale is 
used for display purposes only, ancfthe image has been convolved 
with a (7 = 0.1° Gaussian smoothing kernel. 


We perform a binned maximum likelihood analysi^ 
of the data described above using the following method. 
We bin the gamma-ray photons detected by Fermi into 
0°.l X 0.1° cells, thereby creating a 200 x 200 pixel map 
of the region of interest. Front-converting photons of en¬ 
ergies between 200 MeV and 20 GeV are placed into six 
logarithmically spaced energy bins for fitting of a power- 
law spectrum. By selecting only front converting pho¬ 
tons, we are able to enhance the angular resolution of our 
data. We count photons only when the LMG is within 

See http://fermi.gsfc.nasa.gov/ssc/data/analysis/ 
scitools/binned_likelihood_tutorial .html for details. 
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100° of zenith so that gamma rays originating from Earth 
do not contaminate the data. Finally, we adopt the 
P7REP.CLEAN_V 15 instrument response function as rec¬ 
ommended by the Fermi Science Support Center. 

Our analysis of the gamma rays from the LMC has used 
numerous maps of the stars and ISM. These maps, in¬ 
troduced below, are shown in arbitrary units because the 
normalizations are ultimately left as free parameters in 
our binned maximum likelihood analysis. Our aim here 
is to familiarize the reader with the various spatial distri¬ 
butions of the LMC’s interstellar components. We refer 
the reader to the references given below for further details 
about each map, and we give a more detailed description 
of the physics these maps represent in § [^ The gamma- 
ray map shown in Figure has coarse spatial binning 
compared with measurements at other wavelengths. We 
regrid the maps presented below to match the 200 x 200 
pixel TAN projection of Figure 

To model gas-dominated (hadronic -|- bremsstrahlung) 
gamma-ray emission, we use the maps shown in Figure!^ 
The top image in Figure is a 656.3 nm map of the 
LMC from the Southern Ha Sky Survey Atlas (SHASSA) 
taken using the El Ena no telescope at Cerro Tololo Inter- 
American Observatory ( [Ganstad et al.|2Q^ ). It has been 


extinction co rrected using the 160 /im image of Meixner 
et al. (2006) (se^Appendix [A| for details). The middle 
image m L'igure is an HI 2Fcm line map of the LMC 
using combined data from the Australi a Telescope Coni - 
pact Array and the Parkes Telescope ( Kim et ^^[2003 ). 
Here, we convert intensity to column density by assum¬ 
ing the n eutral hydrogen em issio n to be optically thin . 
Following Abdo et al. (2010a) and Bernard et al. (2008), 
we assume any presence of a dark neutral medium is 
well traced by the HI map. The bottom image is a map 
of molecular hydrogen as esti mated from the NA NTEN 
CO survey data p erformed by Fukui et al. (|2001 ), where 
Yang et al. (2007) have converted from CO to H 2 column 
using Xco = 5.4 x 10^^ H 2 atoms cm“^ (K km 
We use the images presented in Figure to represent 
the ISRF of the LMC. We selected these images based on 
the global spectral energy distribution (SEP) of the LMC 


(llsrael et al.||2010j |Meixner et al.||2013t |Kim et al.||2013|) 

'I'he SEP peaks m the near infrared at /rm and in 

the far infrared at ^200 /im. The top image of Figure 
shows the distribution of light from low- mass stars at 
1.24 /im . These data are part of the 2MASS (Skrutskie et| 
al.|2006). The bottom image of Figure [^traces infrared- 
emittmg dust at 160 /im. This map was obtained using 
the MIPS instrument on the Spitzer Space Telescope as 
part of the Surveyi ng the Agents of a G alaxy’s Evolution 
(SAGE) program (jMeixner et al.||2006|). 

The data described above represent the distributions 
of targets for the collisional processes that ultimately 
lead to gamma-ray production, i.e., the ISM protons and 
ISRF photons with which cosmic-ray protons and elec¬ 
trons collide. Figure]^ shows one of our cosmic-ray flux 
models: star format ion rates at 6.3 Myr (top ) and 12.5 
Myr (bottom) from |Harris & Zaritsky (2009). The au¬ 
thors generated these maps by adopting an initial mass 
function and by fitting a set of isochrones to data they 
collected with the t/, H, V, I filtering system on the 1 m 
Swope Telescope at Las Campanas Observator y as part 
of the Magellanic Cloud Photometric Survey (Zaritsky 
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Fig. 2.— 9° x 9° (90 x 90 pixel) images of LMC components 
used to model hadronic and bremsstrahlung gamma-ray emission. 
Images are normalized so that their maximum pixel values equal 
one, and the color scales are square root. The box es indicate the 
star-f orming region 30 Dora dus. Toy: Ha em ission (|Gaustad et al.| 
|2001| ). Middle: HI er nission (Kirn et al. j 2003l>. Bottom: ti 9 c olumn 


aensit^ 

jet al.| (| 

r estin 


et al. 

1997| 


oR: 


we present a map of the non-thermal 
component oT^ the 1.4 GHz radio continuum emission in 
the LMG, which traces the distribution of synchrotron- 
emitting cosmic-ray leptons. The original r adio map used 
to gen erate this image was presented in [Hughes et al.| 
(2007) and is a combination of interferometric (ATCA) 
and single-dish (Parkes) data. We estimate and subtract 
the thermal contributi on from the original radio data us¬ 
ing the prescription of [Tabatabaei et al. (2007). We use 
a scaled version of the de-reddened tia rnap presented in 
Figure as an estimate for the thermal component. For 
full details about the construction of the map in Figure 
and for a full resolution image, see Appendix [A} 
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Fig. 3.— 9° x 9° (90 x 90 pixel) images of LMC components used 
to model the interstellar radiation field. The boxes indicate the 
star -forming region 30 Do radus. Top: 2MASS 1.24 p,m J band im¬ 
age (iSkruts kie et al.|2006|>. Bottom: MIPS 160 /rm image (|Meixner| 
let al |2006t 7 -- - 


Fig. 4. — 9° x 9° (90 x 90 pixel) star formation rate maps used for 
our concentrated (see §1^ cosm ic-ray distribution models . Top: 6.3 
Myr star formation rare map ( |Harris fc Zarits^|2009|>. B ottom: 
12.5 Myr star formation rate map ( |llarris ^ Ziaritsky|2009| >. 


3. BACKGROUNDS AND FOREGROUNDS 


The unprocessed data shown in Figure include not 
only diffuse emission from the LMC but also Galac¬ 
tic and extragalactic diffuse emission as well as back¬ 
ground gamma-ray point sources. These contaminat¬ 
ing sources must be removed in order to analyze the 
diffuse gamma-ray emission from the LMC. We use 
the iso_clean_front_v05.txt model supplied with the 
Fermi Science Tools package for the isotropic spectral 
template in our analysis. Similarly, Fermi supplies the 
Milky Way model as gll_iem_v05_revl. f its. For each 
of these models, the likelihood estimator fits a multi¬ 
plicative normalization factor that scales the spectrum 
defined by these two models. We also include the four¬ 
teen point sources in our region of interest reported by 
the Fermi Lyge A rea Telescope Second Source Catalog 
(iNolan et al.|2012|), as well as one SIMBAD blazar, PKS 
JU352-6831, which appears to have flared since the re¬ 


lease of 2FGL. Most point sources are modeled using a 
power-law spectrum with two parameters: a normaliza¬ 
tion that scales the spectrum and a power law spectral 
index. The normalization factors of these point sources 
depend on their apparent brightnesses, and the spectral 
indices range from —1.3 to —2.9. We model the bright¬ 
est point source, 2FGL J0601-7037, using a log parabola 
spectrum parameterized by a normalization and two in¬ 
dices describing the exponential energy dependence. The 
top image of Figure shows our background model as fit 
by the binned likelihood analysis, and the bottom image 
shows the background subtracted 9° x 9° central region. 



0.0 0 .2 0.4 0.6 0.8 1.0 

y /counts /max(counts) 

Fig. 5.— Non-thermal 1.4 GHz map of the LMG used for our 
diffuse cosmic-ray electron models. See |Hughes et al.| (|2007|) for 
details about the original radio data. For tuii resolution and mfor- 
mation about the subtraction of the thermal component is given 
in Appendix!^ 

Thin, solid contours indicate the quantity 


\ni - mil 
^/m 


(4) 


where is the number of counts predicted by the back¬ 
ground model in the ith pixel, and rii is the number of 
counts observed in that pixel. These contours are placed 
at 1, 2, 3, 4, and 5cr. The la contour indicates the region 
used for our time variability analysis described below. 
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RA 02000) 

Fig. 6.— Top: LMC background model including fifteen point 
sources, Galactic foreground emission, and isotropic, diffuse emis¬ 
sion. This model was fit to the data in Figure 1^ using a 
binned likelihood analysis with the Fermi Science Tools soft¬ 
ware. Dashed box is the same as in Figure ^ and square root 
color scale is for display purposes only. Bottom: background 
subtracted counts map. The color bar indicates the quantity 
sign(data — background)-^!data — background!. The thin, solid 
contours show cr = |data — background!/^background and are 
placed at 1, 2, 3, 4, and 5cr. The extent of this image is the same 
as the dashed box in the top image. 


4. TIME VARIABILITY 

If the gamma-ray emission from the LMC is indeed 
caused by cosmic-ray energy losses, the gamma-ray flux 
should not vary app reciably on the time scales probed 
by Fermi. Following Abdo et ah (2010a), we perform a 
month-by-month analysis of our data sample to ensure 
the LMC data are not significantly contaminated, for 
example, by a background blazar. 

In the ir analysis of the first year LMC data, 

(2010a) found an unusually high gamma-ray f 


Abdo et al. 
ux for pho- 

tons >100 MeV during the fourth month of observations 
and excluded this time period from their analysis. We 
also find enhanced flux in month four of our time series 
analysis of 200 MeV - 20 GeV photons. Although our 
method of estimating flux differs from that of the Fermi 
study, we also agree that the enhanced emission origi¬ 
nates in or near the 30 Doradus (30 Dor) star-forming 
region. 

Because of the strict requirements we have used for 
event selection (front-converting, zenith angle <100°, 
etc.), we ha ve fewer photons pe r month in our data sam¬ 
ple than didjAbdo et al.|(|2010a|). To avoid dependency on 
a spatial and spectral model tor estimating monthly LMC 


fluxes, we choose to perform photon counting within the 
Icr contour shown in Figure Photon counts are con¬ 
verted to fluxes by dividing the monthly exposures calcu¬ 
lated using the gtexpmap2 function of the Fermi Science 
Tools software. It is important to note that the fluxes 
estimated from this method are sums of the LMC, diffuse 
background, and Milky Way components. This does not 
affect our search for variability because, as for the LMC, 
these diffuse components will also not vary appreciably 
on the time scale of a month. 

Figure shows the results of our time series analysis. 
Month four is the only month with flux that lies out¬ 
side three standard deviations of the mean. Month four 
comprises ~ 1/66 of our data, and we find that removing 
this month from our data would reduce the mean flux by 
0.7%. Compare this to the Fermi study in which month 
four accounted for 1/11 of their data, and by our cal¬ 
culation, affected their flux estimate by 4%. Given its 
< 1% effect, and to avoid complicating the data reduc¬ 
tion process, we choose to keep month four in our data 
set. 

We repeat our time series analysis using the 2cr contour 
shown in Figure which encompasses the 30 Dor region 
and extends to the northwest. Again, the month four 
flux is a unique outlier greater than three standard de¬ 
viations from the mean mon thly flux. This supports the 
claim by Abdo et al. (2010a) that the enhanced emission 
originated in or near 30 Dor. 


6.5 

to 



^ ^'^0 10 20 30 40 50 60 


month 

Fig. 7. — LMC light curve determined by counting photons within 
the thick, dashed contour in Figurej^ The solid black line indicates 
the average monthly flux, the dasned lines indicate Icr, and the 
dotted lines indicate Scr. Month four is the only month with flux 
outside 3cr from the mean. Note: these flux estimates include 
foreground and background photons and are, therefore, higher than 
fluxes reported from our likelihood analysis. 


5. SPATIAL ANALYSIS: MODELS 

Our strategy is (i) to identify physically motivated 
models for the spatial character of the LMG gamma- 
ray emission, (ii) to connect these models to LMG ob¬ 
servables at non-gamma-ray wavelengths, and (hi) to use 
spatial fits of the Fermi LMG signal to test these mod¬ 
els. We perform a morphological analysis of the Fermi 
gamma-ray map of the LMG using several combinations 
of the images in Figures ii and These model 
images are then included in our binned likelihood analy¬ 
sis (described in § |^ as sources of type “DiffuseSource” 
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and spatial models of type “SpatialMap.” The normal¬ 
ization and the single-power-law spectral index are left 
as free parameters of each model map. Our models take 
the functional form 

I^{x,y,E) = '^aiMi{x,y)E~'’% (5) 


where and bi are the free parameters of the fit, Mi 
is a combination of model maps described below, and 
is the LMC gamma-ray intensity as a function of 
plane-of-sky position (x, y) and energy E. The power- 
law dependence of intensity on energy well-describes the 
gamma-ray spectrum at >1 GeV, but at lower energies, 
the spectrum turns over (see ^). The index given by 
Equation ^ is an average over the entire energy range 
and will tend to underestimate the high energy power 
law behavior. 

The models we choose from the maps in Fig nres [^[^[^ 
an d[5]fall into thre e categories: gas models (§ |5.2| ), radia¬ 
tion models (§|5.3|), and a gas + radiation model (§ |5.4[ ). 
Additionally, we subclassify these models as diffuse and 
concentrated based on the quantity used to represent the 
cosmic-ray flux (further explanation below). 

We perform our analysis with and without data from 
the highly active star-forming region 30 Doradus. The 
30 Dor region dominates the gamma-ray, Ha, 6.3 Myr 
star formation, and synchrotron maps, and we remove it 
by subtracting all counts in the 1° x 1° region from the 
model maps. We then model the 30 Dor gamma-rays 
separately as a point source. The region we remove is 
indicated by the boxes in Figures and[^ 


5.1. Gamma-ray Intensity 

The high-energy emission processes we consider are 
from the interaction i + j 7 + • • • of cosmic-ray 
particles i G (e“, ions) with interstellar target j G 
(gas, radiation). The emissivity from such processes is 
schematically given by 


n (F ) = 

dVdtdE. 




( 6 ) 


where is the flux of cosmic-ray species i, and Uj is the 
number density of target species j. For a rigorous d eriva- 
tion of the gamma-ray emissivity, see Stecker (1971). The 
cross section a for gamma-ray production depends on the 
cosmic-ray energy and on the ambient photon energy in 
the case of inverse Compton. Thus, the emissivity aver¬ 
ages over the appropriate cosmic-ray spectrum and ISRF 
as indicated by angle brackets. Our spatial analysis does 
not integrate over energy dependence and instead uses 
observable proxies for cosmic-ray flux and target densi¬ 
ties. 

Because the LMC i s optically th i n to gamma rays 
at Fe rmi energies (e.g. Stecker 1971 Moskalenko et al. 
2006), the photon number intensity or surface brightness 
resulting from emission process I is just 


h{E^) = I qe{E^) dz, (7) 

^los 

the line of sight integral of the emissivity, where we take 
the 2 : axis to be along the sightline. Substituting our 
expression for emissivity, the intensity is a line integral 


of the product of cosmic-ray flux and target density: 

Ii= {^iUjaij^i) dz 

Jlos 

= I dNj = Nj . (8) 

Jlos 

We see that the intensity amounts to a weighted av¬ 
erage of cosmic-ray flux times target column density 
Nj = f Uj dz. Thus, we expect the gamma emissiv¬ 
ity to trace the column density of relevant targets. In 
the case of gas targets, this is the column density of the 
gas phases in which the cosmic rays reside. In the case of 
radiation targets (inverse Compton), the target photon 
column density is, for an optically thin emission region, 
proportional to the observed surface brightness. 

From Equation ([^, we can define two limiting cases 
based on the cosmic-ray flux distribution In one 
limit, does not vary spatially, which amounts to as¬ 
suming diffusion produces homogeneity in the cosmic 
rays and erases memory of source positions. In this case, 
the target distributions uniquely trace the gamma-ray 
distribution. In the other limit, diffusion has not yet oc¬ 
curred, and is traced by sites of cosmic-ray sourcing. 
We describe these limits in the context of our spatial 
models below. 

If the cosmic-ray spectra vary spatially—as is surely 
the case at some level—then the gamma-ray spectrum 
will vary spatially for each emission process. However, 
we will ignore this effect both for simplicity purposes and 
because of the finite amount of spatial information and 
spectral information in the Fermi data. Previo us studies 
of the LMC have made the sa me assumption (Abdo et 
aL]|2QlQa Murphy et al.||2Q12 ). 


5.2. Gas Models 

In the case of hadronic and bremsstrahlung gamma-ray 
emission, the target particles are interstellar gas atoms 
and molecules. The gamma-ray signal from these inter¬ 
actions traces the overlap of cosmic-ray ions and elec¬ 
trons with the LMC gas. In the limit where cosmic- 
ray flux is spatially uniform, we expect the hadronic and 
bremsstrahlung gamma-ray spatial signal to trace the to¬ 
tal gas density. On the other hand, in the limit where 
cosmic rays are confined near their sources, we expect a 
gamma-ray signal concentrated near sites of recent star 
formation. 

Figure shows maps of ionized, neutral, and molecu¬ 
lar hydrogen. These three maps do not trace the thermal 
proton density in the same way. The HI and H 2 maps are 

proportional to column density, A^hi = /q nnds ^ uhF, 
where L is the path length along the line of site. The Ha 
emission, on the other hand, results from recombination 
and is proportional to UpUeds which is essen¬ 

tially density-weighted column density. Since nn ranges 
from below 0.1 cm“^ in the diffuse medium to over 10^ 
cm“^, using the Ha surface brightness as a proxy for the 
ionized hydrogen column density can produce differential 
errors by factors greater than 10^. We suggest the use 
of the square root of Ha surface brightness as a proxy of 
ionized column density. The square root of Ha emission 
is proportional to Because the emission path 

length ranges from 10 pc for individual HII region to 500 
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pc for the diffuse ionized gas (DIG) in the disk, the differ¬ 
ential error introduced by our under-counting the path 
length is less than a factor of 10. Thus, while not ideal, 
we consider the square root of Ra as a far better proxy 
than Ra for the ionized gas column density. 

Our first gas models use the HI, H 2 , Ra (included 
for comparison purposes) and maps by themselves. 
These mode l s, the first three of which are used in the 
Abdo et al.| (|2010a|) study, are akin to assuming an ex¬ 
treme case in which cosmic rays have diffusively propa¬ 
gated away from their acceleration sites in such a way 
as to create a uniform flux distribution across the LMC. 
With no spatial dependence, the isotropic cosmic-ray flux 
level then becomes one factor in the normalization pa¬ 
rameter, a^, in Equation <§■ 

We explore the opposite extreme, a case in which cos¬ 
mic rays remain in the immediate vicinity of their accel¬ 
eration sites, by using the maps of star formation rates in 
Figure as proxies for the cosmic-ray flux. These mod¬ 
els are constructed such that the hydrogen maps (HI, 
H 2 , and \/H(a) are multiplied by each of the star forma¬ 
tion rate maps, so that the binned likelihood software is 
given two separate functions with which to fit the LMC 
gamma rays. For example, our Hl-times-star formation- 
rate model is formed by multiplying the HI map by both 
the 6.3 and 12.5 Myr star formation rate maps, i.e., 

Mi{x,y) = {/m'^e.SMyr, 7HI'012.5Myr} , (9) 

where V^e-SMyr and ' 0 i 2 . 5 Myr are the two epochs of star 
formation. These two images are then simultaneously fit 
as separate components, each with its own normalization 
and spectral index. This method allows us to comment 
on the relative importance of the 6.3 and 12.5 Myr star 
formation epochs by comparing their relative flux contri¬ 
butions to each model. We combine the two maps fit on 
this first pass into one by adding them based on their flux 
contributions. After this step, the concentrated models 
take the form 

7 \/r/ \ -^6.3Myr j , -^12.5Myr j , 

^ -7HI'06.3Myr H- ^ -7HlV^12.5Myr• 


^tot 


tot 


( 10 ) 

Finally, this single map is then refit to the data using 
a single normalization and a single spectral index pair. 
Maintaining the same number of free parameters for both 
diffuse and concentrated models eases their comparison 
using the likelihood statistic described below. 

Conclusions about models that fit multiple images each 
with its own normalizations must be made cautiously; 
any physical offset between the images is washed out by 
the free renormalization p erformed by the like lihood esti¬ 
mator. For example, when Abdo et al. (2010a), construct 
their A/'(HI) + N{R 2 ) + A/'(HH) model where they allow 
each hydrogen component to be fit independently, some 
information is lost. Specifically, the fact that HI typically 
has the greatest column density along a given site line is 
disregarded. Fortunately, the total LMC star formation 
rate is ^0.4 M^:^ yr~^ during bot h the 6.3 and 12.5 Myr 


epochs (Harris & Zaritsky 2009, Figure 11), so the two 
maps have virtually no physical offset. 

We exclude from our discussion the older star for- 


might suggest that cosmic rays accelerated during that 
epoch have undergone significant diffusion. The LMC 
was forming stars nearly 25% more slowly 25 Myr ago 
( Harris fc Zaritsky]|2009| Figure 11). 

hbr completeness, we have also constructed a model of 
the form 

Afi=/HI+/H2 5 (11) 

as a “total” hydrogen-mass model (note this model is 
allowed only one normalization and one spectral index 
parameter). We have excluded Ra from this model, first 
because /hq does not measure the same quantity as /hi 
and /h 2 and also because ionized hydrogen does not con¬ 
tribute s ignificantly to the total hydrogen mass (Yang et 
al.|2007D . 


mation rate maps of Harris & Zaritsky (2009) because 
when we included the 25 IVlyr map, it was not given an 
appreciable weighting by the likelihood analysis. This 


5.3. Radiation Models 

Inverse Compton emission is the result of scatterings 
between cosmic-ray electrons and photons of the ISRF. 
As with the gas models, we split the radiation models into 
diffuse and concentrated categories. Again, we use the 
star formation rate maps as a measure of the cosmic-ray 
flux in the concentrated models. We use the map pre- 
sented in Figure fo r our diffuse models. The ^3 GeV 
(Murphy et al.|20^ ) electrons emitting synchrotron radi- 
ation at 1.4 GHz are not themselves energetic enough to 
significantly contribute to the inverse Compton gamma- 
ray emission. However, this population of lower-energy 
electrons will be partially comprised of formerly high- 
energy electrons that have undergone inverse Compton 
scattering. These older, lower-energy electrons have had 
additional time to propagate after their scattering events 
and should be more diffuse. 

Our choices in models for the ISRF are mot i vated by 
local maxim a in t he SEDs of Israel et al. (2010), Meixner 
et al. (12013), and Kim et al. (2013|). We use the 1.24 /am 
map shown in Figure which traces the dominant en¬ 
ergy contribution from stars, and we use the 160 /im 
map as a tracer for the energy contribution from dust. 
Although the cosmic microwave background (CMB) is 
not included in the SEDs mentioned above, it neverthe¬ 
less contributes a relevant fraction of energy to the ISRF. 
The CMB is an isotropic source, so for our CMB models, 
in practice, we fit either the 1.4 GHz map or the star 
formation maps by themselves. 

We note that the minimum electron energy required 
to scatter a 1.24 ym photon up to 200 MeV is given 

by Ee = 7minTOeC^ = 0.5ye^/ei.24 = 3.6 GeV. 

This is the lower limit energy for a cosmic-ray electron to 
contribute to the Fermi data we analyze here given the 
ISRF energies trace by the maps described above. So our 
previous statement that 1.4 GHz synchrotron-emitting 
cosmic-ray electrons do not contribute significantly to 
inverse Compton is well justified. 

The map of synchrotron radiation traces the magnetic 
field strength in conjunction with the cosmic-ray electron 
density. We have made the simplifying assumption that 
at scales probed by Fermi^ the magnetic field strength 
does not vary appreciably and is therefore a constant 
scale factor across the entire map. 

5.4. Gas + Radiation Model 

We expect both the ISRF and the ISM to contribute 
as targets, and we explore models in which gas and ra- 








































diation maps are fit simultaneously. These models take 
the form 





f Fe.SMyr 
\ Ff 

\ tot,gas 

f Fq_SM yr 

\ -^tot,rad 


V^e.SMyr + 
V^e.SMyr + 


-^12.5Myr 

-^totjgas 

-^12.5Myr 

-^tot,rad 


' 012 . 


5Myr 


012.5Myr 



( 12 ) 


The radiation and gas components are each fit using their 
own normalizations, a^, and spectral indices, 6^, because 
we do not assume that inverse Compton and the collision 
processes involving gas targets result in similarly shaped 
gamma-ray spectra. 

We do not test every combination of the gas and radi¬ 
ation maps described above because of resource limita¬ 
tions. We explain below which maps we have chosen for 
these combined models after discussing the performances 
of all single map models. 


6. SPATIAL ANALYSIS: RESULTS 


As did Ackermann et al. (2012) in the Fermi group’s 
study of the IVlilky Way gamma-ray emission, we use 
the likelihood statistic to distinguish between the spa¬ 
tial models described in the previous section. The likeli¬ 
hood as computed by the Fermi Science Tools softwar^ 
is given as 


LH = e 


-Ne: 


n 


m- 


(13) 


across all of our models are as follows: (1) concentrated 
models, i.e., models that use star formation rate maps 
to represent the cosmic-ray flux, better fit the gamma- 
ray data than do diffuse models a nd lead to conc l usions 
con sistent with the conc lusions of Murphy et al. (2012) 
and Abdo et al. (2010a|). (2) Spectral indices ht to our 
subset of the LiVlC Fermi observations consistently lie 
between 2.2 and 2.4 with the exception of the 160 ym 
diffuse model, which we address below. This is the case 
whether 30 Dor is included or not. 


The total LMC flux >200 MeV estimated from our 
models range from 0.92 x 10“^ to 1.96 x 10“^ ph cm“^ s“^ 
for the cases in which we include the 30 Doradus data, 
and fluxes range from 0.94 x 10“^ to 1.40 x 10“^ ph cm“^ 
s~^ when we exclude 30 Doradus. As expected, for each 
given model, estimated fluxes are higher when we include 
the 30 Doradus data, with the one anomaly being the 
diffuse 160 /im radiation model. Spectral indices range 
from 2.22 to 2.52. 

Model-to-model changes in the LMC flux values are 
mirrored by compensating changes in the estimated 
Milky Way flux. The difference between the LMC and 
the Milky Way in the absolute variation in flux is ex¬ 
plained by the fact that the Milky Way’s flux is across 
the entire sky, whereas our region of interest is ^1% of 
the sky. The flux differences between models with 30 Dor 
and those without are well explained by the 30 Dor point 
source flux in conjunction with the Milky Way flux. 


where i is an index over image pixels in both space and 
energy, indicates the number of counts predicted by 
the model at pixel i, is the observed number of counts 
at pixel 0 and Aexp is the total number of observed 
counts. We rank the spatial models based on their asso¬ 
ciated likelihood statistic. 

To quantify the robustness of our ranking method, we 
generate 10 boots trapped realizations of our data (e.g. 
Press et al. 2007). Given the N photons in the Fermi 
data set, our method is to select random samples of size 
N photons allowing for replacement, i.e., some individual 
photons from the original data set are used more than 
once in a given realization, others not at all. We have 
released our code to generate bootstrapped samples of 
the Fermi data at http://github.com/garyForeman/ 
FermiBootstrap, 

We hnd a wide spread in LH of the null hypothesis 
(i.e. the background model) among our bootstrap sam¬ 
ples, and therefore we are not able to comment on the 
significance of the the likelihood statistic as a means of 
ranking LMC spatial models. What we can report is the 
ranking of each model as used for the true data set and 
the frequency with which the same ranking is reproduced 
among the bootstrapped samples. In Figuresandthe 
error bars in the upper left plot indicate the mean and la 
spread of the ranks across our 10 bootstrapped samples. 
We will explore more advanced statistical techniques in 
future work. 

We present the results of our binned likelihood anal¬ 
yses in Tables mi The broad statements we can make 

^ See http://http://fermi.gsfc.nasa.gov/ssc/data/ 
analysis/documentation/Cicerone/Cicerone_Likelihood/ 
Likelihood_formula.html for details. 


6.1. Gas Models 


As did I Abdo et al. 


(2010a), we find the Ha model 
consistently outranks the HI and H 2 models both when 
including and excluding the 30 Dor region. However, 
we find that that our optically thin HI model better fits 
the gamma-ray data than H 2 whereas the Fermi study 
found the inverse. We are unfortunately preclu ded from 
ma king a m ore quantitative comparison because [Abdo et| 


al. (2010a) do not report the normalization and spectral 
index parameters estimated from their maximum likeli¬ 
hood fits of these spatial templates. Our stricter require¬ 
ments on photon quality and the extra point sources we 
have used in our background model exacerbate this dif¬ 
ficulty. 

The upper-left plot of Figure reveals that when 30 
Dor is included, concentrated models are consistently 
ranked higher than their diffuse counterparts. We remind 
the reader that diffuse gas models assume a spatially uni¬ 
form cosmic-ray flux distribution, whereas concentrated 
models approximate the cosmic-ray flux using the star 
formation rate maps of Figure The four best-fitting 
models are indistinguishable based on our bootstrapped 
ranking errors, but, again, we do not consider the diffuse 
Ha model to be a good estimate of the ionized hydrogen 
column density. This model has bee n included only for 
comparison with Abdo et al. (2010a). 

When we exclude the 30 Dor region, we find the diffuse 
Ha and \/Ha maps to consistently outrank all concen¬ 
trated models. We remind the reader that we do not 
consider the Ha model to be physical because Ha inten¬ 
sity is not proportional to the ionized hydrogen column 
density. However, the high ranking of the VHa does 
stand out, especially when comparing its In(LH) value 
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TABLE 1 
Gas Models 



LMG 

Model 

ln(LH)+ 

1.2 X 10® 

LMG Flux 

200 MeV - 20 GeV 
(lO”’^ ph cm“^ s“^) 

Spectral 

Index 

ai 

a2 


H2 

-8250.4 

1.58 ±0.03 

2.32 ±0.02 




HI 

-8176.8 

1.96 ±0.03 

2.27 ±0.02 




HI+H2 

-8107.2 

1.96 ±0.3 

2.27 ±0.02 



With 

H2 (ai 'IpG.S + CL2 '012.5) 

-7745.8 

1.48 ± 0.02 

2.34 ±0.02 

0.663 

0.337 

30 Dor 

Vh^ 

-7569.8 

I.9I±0.03 

2.22 ±0.01 




Ha 

-7413.7 

1.53 ±0.02 

2.29 ±0.02 




\/Ha (ai 06.3 + CL2 012.5) 

-7388.7 

1.49 ± 0.02 

2.29 ±0.02 

0.204 

0.796 


hi (ai 06.3 + a2 012.5) 

-7388.6 

1.57 ±0.02 

2.28 ±0.02 

0.193 

0.807 


HI+H2 (ui 06.3 + «2 '012.5) 

-7378.0 

1.57 ±0.02 

2.28 ±0.02 

0.207 

0.793 


H2 

-7750.8 

I.II±0.03 

2.34 ±0.02 




HI 

-7700.6 

1.40 ±0.03 

2.28 ±0.02 




HI+H2 

-7673.9 

I.4I±0.03 

2.27 ±0.02 



Without 

H2 {ai 06.3 + 0,2 012.5) 

-7631.3 

1.02 ±0.03 

2.34 ±0.02 

0.569 

0.431 

30 Dor 

Vh^ 

-7383.7 

1.40 ±0.03 

2.22 ±0.02 




Ha 

-7336.9 

1.25 ±0.03 

2.26 ±0.02 




VHo {ai 06.3 + Ci2 012.5) 

-7532.5 

1.02 ±0.03 

2.34 ±0.02 

1.000 

0.000 


HI (ai 06.3 + CL2 012.5) 

-7409.3 

I.I4±0.03 

2.26 ± 0.04 

0.209 

0.791 


HI+H2 {ai 06.3 + «2 012.5) 

-7395.6 

I.I4±0.03 

2.26 ±0.02 

0.205 

0.795 

Note. — 

Models ordered by their with 30 Dor 

ranking based on In(LH). Here, ai and 02 represent the fraction of flux contributed by the 6.3 


and 12.5 Myr star formation rate maps respectively. 


TABLE 2 

Radiation Models 



LMG 

Model 

ln(LH)+ 

1.2 X 

LMG Flux 

200 MeV - 20 GeV 
(10“^ ph cm“^ s“^) 

Spectral 

Index 

ai 

a2 


1.4 GHz X 160 iim 

-8651.4 

0.92 ± 0.02 

2.52 ± 0.03 




1.4 GHz X 1.24 /im 

-7862.4 

1.34 ±0.02 

2.29 ± 0.02 



With 

1.24 /im(ai 06.3 ± cl2 012.5) 

-7569.7 

1.46 ±0.02 

2.31 ±0.02 

0.977 

0.033 

30 Dor 

GMB (ai 06.3 ± a2 0i2.5) 

-7549.1 

1.63 ±0.03 

2.26 ± 0.02 

0.839 

0.161 


1.4 GHz X GMB 

-7463.1 

1.60 ±0.03 

2.26 ± 0.02 




160 iiicn{ai 06.3 ± cl2 012.5) 

-7408.0 

1.51 ±0.02 

2.28 ± 0.02 

0.174 

0.826 


1.4 GHz X 160 fim 

-7595.2 

1.23 ±0.03 

2.34 ±0.02 




1.4 GHz X 1.24 /im 

-7895.6 

0.94 ±0.03 

2.29 ± 0.03 



1 1" In 011 1 

1.24 I2m(ai 06.3 ± ^2 012.5) 

-7513.1 

1.05 ±0.03 

2.31 ±0.02 

1.00 

0.00 

V V 1 UiilJ u. t 

30 Dor 

GMB (ai 06.3 ± CL2 012.5) 

-7396.0 

1.11 ±0.03 

2.25 ± 0.02 

0.391 

0.609 

1.4 GHz X GMB 

-7460.6 

1.34 ±0.03 

2.27 ±0.02 




160 06.3 ± a2 012.5) 

-7611.1 

0.98 ± 0.03 

2.36 ± 0.03 

0.588 

0.412 


Note. — Models ordered by their with 30 Dor ranking based on In(LH). Here, ai and a 2 represent the fraction of flux contributed by the 6.3 
and 12.5 Myr star formation rate maps respectively. 


TABLE 3 

Gas + Radiation Model 



LMG 

Model 


ln(LH)+ 

1.2 X 

LMG Flux 

200 MeV - 20 GeV 
(10“^ ph cm“^ s“^) 

Spectral 

Index 

ai 

a2 

With 

30 Dor 

ai\/Ha ± 

a2 160 )um(C06.3 ± (1 

- C) 012.5) 

-7307.4 

1.72 ±0.08 

2.22 ± 0.05 
2.26 ± 0.03 

0.449 

0.551 

Without 

30 Dor 

aiVaoi ± 

(22 160 /2m(C'06.3 ± (1 

— C) 012.5) 

-7321.9 

1.35 ±0.06 

2.21 ±0.03 
2.28 ±0.07 

0.761 

0.239 


Note. — Likelihood statistic reported here not directly comparable to Tables [T] and because more free parameters are used in these models. 
Here, ai and 02 represent the fraction of flux contributed by the gas and radiation maps respectively. C and (1 — C) correspond to a\ and 02 from 
Table [ 2 ] 
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Fig. 8.— Bootstrapped results for the ^s models. Triangles (full LMC) and squares (30 Dor removed) indicate the values found from 
the Fermi data set (data listed in Table [l]> . The error bars indicate a Tier spread on the mean quantities found from ten bootstrapped 
samples. We have slightly offset the full EMC data points from the 30 Dor removed data points to avoid error bar confusion. Gray shaded 
regions are meant to help distinguish between quantities from separate models. Upper left: rank as determined by the likelihood statistic. 
Upper right: LMC ffux. Lower left: spectral index. Lower right: Milky Way foreground ffux. 
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Fig. 9.— Same as Figurebut for radiation models. Note: the best fitting gas model, “HI+H2, C,” consistently outranks the “160 /rm, 
C” model. 


with the other models (see Table . 

The LMC flux varies by ^25% depending on which 
spatial model is used to to estimate the flux. The high 
fluxes from the diffuse HI and diffuse HI+H2 models can 
be understood by the fact that the the HI emission ex¬ 
tends farther from the LMC center than do emissions 
from other components. Spectral indices vary at about 
5% depending on the spatial model, and the estimated 
Milky Way fluxes vary at about 4%. 

6.2. Radiation Models 

For the radiation models plotted in Figure the state¬ 
ment that concentrated models better explain the LMC 
gamma-ray data is less strong. This is a reflection of the 
different diffuse model we use to represent the cosmic- 
ray electron flux. Rather than assuming a spatially uni¬ 
form distribution, we use the synchrotron radiation map 
presented in Figure ^ The synchrotron map traces low- 
energy cosmic-ray elS:trons, some of which have lost en¬ 
ergy through inverse Compton scattering. These older 
electrons have had more time to diffuse throughout the 
LMC than their younger, high-energy counterparts, so 
we use them for our diffuse model. 


The diffuse 160 /im model exhibits strange behavior in 
that it becomes brighter when the 30 Dor photons are 
removed. This finding is explained by the fact that this 
model is highly peaked in the 30 Dor region, which in 
turn dominates the fit. The rest of the diffuse emission 
is under-explained by this model and is ultimately in¬ 
corporated into the fit of the diffuse Milky Way model 
(see lower right plot of Figure]^. When the peaked 30 
Dor region is removed, the rest of the LMC photons are 
better explained by the model resulting in a higher LMC 
flux estimate. 

The spread in estimated model parameters from the 
radiation models is heavily influenced by the diffuse 160 
jam model, which we consider anomalous. If we include 
the 160 jam model, the spread in LMC flux is about 44%; 
the spread in spectral index is ^10%; and the spread in 
Milky Way flux is ^5%. These spreads reduce to respec¬ 
tively 18%, 5%, and 3% if we discount the diffuse 160 
jam model. 

Finally, we wish to note that the concentrated HI+H 2 
model consistently outranks the best fitting radiation 
model. 
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6.3. Gas + Radiation Models 


As mentioned in §5.4[ we have chosen several combi¬ 
nations of best fitting gas and radiation models to be 
fit simultaneously. Among these models, we find the in¬ 
verse Compton contribution to the fiux to be generally 
between 30% and 50%. We choose one fiducial model to 
present here, which is the diffuse \/Ha model combined 
with the concentrated 160 /am. This model consistently 
outranked other combinations. Results of the model fits 
are listed in Tabled 

This model estimates the inverse Compton flux contri¬ 
bution at 55% when the 30 Dor region is included and 
at 24% when 30 Dor is excluded. The spectral indices fit 
to the gas and radiation components are consistent with 
each other for both data sets. Note the improvement in 
likelihood of the combined model compared with the dif¬ 
fuse VHof model. This likelihood improvement accompa¬ 
nies a lower estimated LMC flux for the gas + radiation 
model. 

In summary, we find that the best fit gas model con¬ 
sistently outranks the best fit radiation model. How¬ 
ever, when we combine gas and radiation models and 
fit simultaneously, we find significant flux is contributed 
by the radiation component, which suggests the leptonic 
contribution to the LMC gamma-r ay si gnal should not 


old for neutral pion production (T’thr = 1-2 GeV). The 
second free parameter is the high-energy power law in¬ 
dex, a. The proton-proton collision cross section is given 
as cTpp = 32(0.96 + g4.4-2.4a^ mbarn, and the exponent 

5 = Q.14(a- ^-^ + 0.44. _ _ 

Similarly, jChakraborty & Fields (2013) give an inverse 
Compton fitting function tor nornial, i.e., non starburst¬ 
ing, galaxies: 


logic ( ) = log 


dt dEJ VV^MW J (15) 

[-40.26 + X(0.0949 + X(0.0503 + 0.0251X))], 

where X = log]^o(^/GeV), and ^ is the star formation 
rate. The functional form of this inverse Compton spec¬ 
trum is based on a one-zone model of the Milky Way that 
assumes electron calorimetry, which leads to linear scal¬ 
ing with star formation rate. Equation is normal¬ 
ized to the total Milky Way inverse Compton lu minosity 
found in th e GALPROP plain diffusion model ([Strong 
et al.||2QlQ ). Because of this normalization. Equation [T^ 
implicitly assumes the same cosmic-ray acceleration eim 
ciency. Other cosmic ray models (e.g., diffusive reaccel¬ 
eration, different confinement volumes) can lead to dif¬ 
ferent efficiencies. Equation also assumes that the 
galaxy with star formation rate has the same ratio of 


be ignored (further discussion in § 8.2). We also find 


concentrated models generally outperform their diffuse electron energy-loss rates (i.e. 5ic/5sync and 5ic/^b] 


c ount erparts, favoring short diffusion scale lengths (see 

§[ 0 ). 


7. SPECTRAL ANALYSIS 

In this section, we consider the spatially integrated 
LMG gamma-ray spectrum from the data presented in 
Eigure We split our data into six logarithmically 
spaced energy bins from 200 MeV to 20 GeV and repeat 
the binned likelihood analysis on each bin independently. 
To estimate the gamma-ray flux from each bin, we use 
the diffuse VHa + concentrated 160 /am model. In this 
analysis, we model the LMG gamma-ray emission with 
twelve free parameters: two normalizations per energy 
bin, where the spectral index within each energy bin has 
been fixed to zero. As for the LMG models, all point 
sources have been fixed to have a spectral index of zero. 

Performing our analysis on individual energy bins 
in this way allows us to remove the assumption of a 
single power law used in our spatial analysis. This 
methodology also allows us to investigate the pion decay, 
bremsstrahlung, and inverse Gompton flux contributions 
by fitt ing theoretical curves to the LMC gamma-ray spec¬ 
trum. jPfrommer & EnBlin (2004) give a fitting function 


as the Milky Way. This is not true of the LMC, which 
has a much higher neutral hydrogen density and a lower 
ISRE energy density. Therefore, we scale Equation (15) 
by an additive factor of 


log 


10 


^ ^LMC,ISRF ^tot,Mw(l GeV) \ 

V ^MW,ISRF ^tot,LMc(l GeV) J 


- 1 , 


(16) 


for the gamma-ray spectrum predicted by pion decay: 




_ = + rT rN ^ 

dAdtdEdn 4n ^ GeV 3a 


GeV 


2E^ 


m^oc^ 


+ 


2E^ 


m^oc^ 


-S' 


oijS 


(14) 

The first free parameter of this function is a normal¬ 
ization factor that is proportional to AhihcRp, the hy¬ 
drogen column density times a proxy for the number 
density of cosmic-ray protons above the energy thresh- 


where Visrf is the energy density of the ISRE, and 5tot 
is the sum of four electron energy loss rates: ^sync, 
^hrem ; cind 6ion, whcrc 6ion Is thc lo ss rate from ioniza¬ 
tion (Ginzburg & Syrovatskii 1964). Eor further details 
about the quantities we have used for the LMC, please 
refer to Appendix We adopt Vmw.tsrf. and 

^tot.Mw fl GeV) from the fiducial Chakraborty & Eields 
( 2Q13| ) model. The scaling of Equation (16) is moti¬ 
vated by the unnumb ered relation that appears be tween 
Eqs. (23) and (24) of [Chakraborty & Eieldsj ([2013 ). 

Because of the calorimetry of inverse Compton emis¬ 
sion, Equation (f^ rescaled by (p!^ gives a zero param¬ 
eter prediction for the LMC inverse Compton luminos¬ 
ity, if we assume a universal cosmic-ray acceleration effi¬ 
ciency. The detailed shape of the inverse Compton spec¬ 
trum de pends on the details of the L MC ISRE, but as 
shown in Chakraborty & Eields[ ([2013 ), in the calorimet¬ 
ric limit, the dependence is mild for a plausible range of 
radiation field components (UV, optical, IR, and CMB). 
Alternatively, if one adopts the observed LMC star for¬ 
mation rate, then the inverse Compton luminosity nor¬ 
malization measures the ratio ee,LMc/ee,MW of cosmic- 
ray electron efficiency in the LMC to that in the Milky 
Way. 

We have formulated our own bremsstrahlung fitting 
function based on the same one-zone model discussed in 
the previous paragraph, which is a modified version of 


the Chakraborty & Eields (2013) model. The detailed 


derivation of this fitting function can be found in Ap- 
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pendixj^ The equation we fit takes the final form 


dN^ 
dt dE 


= 5.93 X 10^^ s 


29 -1 V^LMC 4(a 


X < 


dEe 


'0MW E 




(Jt ln(183)nHcK 


-1 

7 


L 


- 0.8 


- 1.25 _ ^- 1.25 


E^ < El) 


)] Eb< < E„ 


(17) 


This spectrum is based on a broken power-law spectrum 
of the cosmic-ray electron energies E^ with break en¬ 
ergy El) = A GeV and a hard cutoff at ^max = 2 TeV. 
Here a is the fine structure constant, (Jt is the Thomson 
scattering cross section, nn is the total hydrogen den¬ 
sity, and b{Ee) is the sum of electron energy-loss rates 
from bremsstrahlung, inverse Compton scattering, syn¬ 
chrotron, and ionization. 

The reason bremsstrahlung and inverse Compton spec¬ 
tra depend on the star formation ratio is that the cosmic- 
ray electron injection rate is proportional to the super¬ 
nova rate, which is proportional to the star formation 
rate. Both Eqs. ([T^ and Eqs. ( p!7[ ) have been normalized 
to the Milky Way spectra calculated from GALPROP 
(Strong et al. 2010). Therefore, to apply these spectra 
to other galaxies with different cosmic-ray densities, one 
scales by the star formation rate. 

Ideally, we could constrain the free parameters of 
Eqs. (pl^, (p!^, and (pT|) using the maximum likelihood 
estimator implementeTin the Fermi Science Tools soft¬ 
ware. Unfortunately, the normalization is the only free 
parameter allowed for user-defined spectral functions, so 
we would not be able to estimate the high-energy spec¬ 
tral index of Equation ( p^ . We proceed by using least 
squares to fit these spectral profiles to the binned flux 
data. The contribution to the fit from each energy bin is 
weighted by the flux error. 

Eigure 10 shows the results of the fits described above. 
Each plot^ows a different subset of the EMC data. The 
solid lines indicate the results of the simultaneous fits 
of Eqs. (p^, (p!^, and (pT|). The dashed, dotted, and 
dotted-dashed Jmes show me contributions from the pion 
bump, bremsstrahlung, and inverse Compton models re¬ 
spectively. We present the fitted parameters in Table 

The LMC gamma-ray flux from 200 MeV to 20 Gev 
is (1.37 ± 0.02) X 10“^ ph cm“^ s“^, which we estimate 
by integrating the fitted spectral function indicated by 
the solid curve in the top panel of Eigure pR| The error 
bounds are calculated by performing the least squares fit 
on an ensemble of bootstrapped samples. 

We can estimate the flux contributions from each of the 
gamma-ray production channels. Eor the full LMC data 
set, our spectral analysis suggests that 6% of the LMC 
gamma-ray flux is caused by inverse Compton scatter¬ 
ing, and 44% is from bremsstrahlung. Eor regions out¬ 
side of 30 Dor, our spectral fit estimates that 3% of the 
gamma-ray flux is from inverse Compton and 18% is from 
bremsstrahlung. Eor the 30 Dor point source model, the 
spectral analysis results in an inverse Compton flux con¬ 
tribution of 12% and a bremsstrahlung contribution of 
80%. 

As indicated in Table the A^nhcRp and the 
'0LMc/'0MW parameters are not particularly well con¬ 
strained. The reason the errors on the spectral flux es¬ 





Fig. 10.— Spectral fits to the bi nned flux estimates. The fit is 
performed using Equations ( |14| > + ( |15| > and three free parameters 
are fit simultaneously: NnhcRp, o;, and bLMc/'bMW- Solid lines 
show the result of the fits, and the dashed, dotted, and dot-dashed 
lines show the contributions from the pion bump, bremsstrahlung, 
and inverse Compton models respectively. Top: model for the full 
LMC. Center: model for the LMC data excluding 30 Dor. Bottom: 
point source model for 30 Dor. 


timate are so small in spite large variation in the fit¬ 
ting parameters is that a Pearson correlation coefficient 
of r = —0.998 indicates NnhcRp and '^lmc/'^mw are 
nearly perfectly anticorrelated. This is caused by the fact 
that the sum of the bremsstrahlung and inverse Comp¬ 
ton spectra as parameterized by V^lmc /'0mw has such a 
similar shape to that of the pion spectrum parameterized 
by A^HhcRp. 

The high-energy power-law index is better constrained 
than the two normalization parameters and consistent 
across all three data sets. An index a = 2.4 is consis¬ 
tent with the result of all three fits. This finding has 































14 


TABLE 4 
Spectral Fits 


Data Set 

A^H^CRp 

max(A^HRCRp) 

OL 

'bLMc/'bMW 

max('0LMc/'bMw) 

Pion Flux 

Brem Flux 

IC Flux 

Full LMC 

1.3 

3.4 

2.4 ±0.3 

0.3 

0.7 

6.9 

6.1 

0.9 

No 30 Dor 

2.1 

2.8 

2.5 ±0.2 

O.I 

0.6 

9.1 

2.0 

0.3 

30 Dor Point Source 

0.03 

0.75 

2.2 ± 1.9 

O.I 

0.1 

0.2 

2.0 

0.3 


Note. — NnhcKp given in units of 10 ^^ cm ^. Fluxes given in units of 10 ® ph cm ^ s ^. max(AHncRp) computed by setting ' 0 lmc/' 0 mw = 0 - 
max(' 0 LMc/' 0 Mw) computed by setting A^HncRp = 0 . 


a strong implication for the underlying cosmic-ray pro¬ 
ton spectrum itself as we discuss below. As mentioned 
above, the spectral indices found from our spatial fits, 
which have a typical value of about 2.3, slightly under¬ 
estimate the high energy power law index a because the 
the spatial models assumed a power law over the entire 
energy range. 

8. DISCUSSION 

8 .1. Cosmic-ray Density and Magnetic Field Strength 

At high energies, the tt^ spectrum (and by exten¬ 
sion the gamma-ray spectrum) is a power law with in¬ 
dex equal t o that of the cosmic-ray protons that creat e 


dex of a = 2.0, which was found by Abdo et al 


them (e.g., Stecker I P^ |d ( 


1986 Gaisser 1991) 


This means the a parameter in Equation! 14) is itself 


the cosmic-ray proton spectral index. By fitting Equa¬ 
tion (14) to the energy-binned gamma-ray data of the 
entire LMC, we find hcRpA^H = 1.3 x 10^^ cm“^, and 
the power- law index of the cos mic-ray proton population 
is a = 2.4. Abdo et al. (2010a) take J N^dQ = 3.6 x 10^^ 
^ - 


H-atom cm ~^ sr a s an average of measurem ents made by 
(2003) and Bernard et al.| ( [2008 ). This result 
corresponds to Nn = 1.4 x 10^^ H-atom cm“^, assuming 
an angular size of 81 deg^ for the LMC (size of Figures S 
1^ and[^. From our fit, hcRn = 9.3 x 10~ ^^ cm“^. 


From ncRp and Cf, Pfrommer & EnBlin (2003) give the 
cosmic-ray proton kinetic energy density as 


^CRp = 


^CRp^pC" / rripC^ 

2{a-l) VGeV 
^ f a — 2 3 — cf 


1—a 


2p 


~1—a 


■p2-l 


(18) 


Here, Bx{a,b) is the incomplete beta function, x = 
(1 and p = Pmin/(^pc). We take Pmin to be 

the minimum cosmic-ray proton momentum necessary 
for pion production, which is 780 MeV c“^ calculated 
from relativistic proton collision kinematics. From Equa¬ 
tion (18), we find the cosmic-ray proton energy density 
to be 3.1 X 10“^^ erg cm“^. If we assume the energy 
density of cosmic rays is dominated by protons and that 
magnetic fields and cosmic rays are in energy equiparti- 
tion, the result corresponds to a magnetic field strength 
of 2.8 pG. 

Our mag netic field strength esti mate is in good agree¬ 
ment with Gaensler et al. (2005), who found the field 
strength to be 4.3 /iG based on Laraday rotation mea¬ 
surements of the LMC. Mao et al. (2012) also calculated 
an equipartition-based ma gnet icneld s t rength of ^2 pG 


based on th e results fro m Abdo et al. (2010a). In their 


calculation, Mao et al. (2012 | ) used a harder spectral in- 


(2010a) 


by fitting a power law plus exponential cutoh to the 
gamma-ray spectrum. This means the power-law index 
they report fits the low-energy end of the spectrum, the 
regime of the pion bump turnover. As mentioned above, 
the high-energy portion of the gamma-ray spectrum is 
determined by the cos mic-ray proton spect r um, w hich is 
better modeled by the [Pfrommer fc Enfilin (2004) fitting 
function. Equation ( p^ 

Despite the discrepancy in details, we are able to draw 
th e same conclusion s from our magnetic field calculation 
(2012). Those authors made a second esti- 


Mao et al. 


mate of ^7 /iG using radio synchrotron data. The agree¬ 
ment between the equipartition and synchro tron esti¬ 
mates along wi th the measureme nt reported by jGaenslCT 


et al. (2005) led Mao et al. (2012) to conclude the cosmic- 
ray energy equipartition with niagnetic fields is not vio¬ 
lated in the LMG. 

8 .2. Leptonic Gamma-ray Emission in the LMC 


The 

previous 

et al. 

1 

2010 a 



revious studies of the the Fer mi LMG data (Abdo 
Murphy et al. 2012) assumed that the 


gamma-ray eniission is dominated by hadronic collisions 
and subsequent neutral pion decay. Both our spatial and 
spectral analyses suggest leptonic processes contribute 
significantly to the LMC gamma-ray emission. 

Our spectral analysis suggests the leptonic contribu¬ 
tion to the spatially integrated LMC spectrum is at 
^50%. Admittedly, the parameters fit in our spectral 
analysis are not well constrained, but the estimated value 
of the LMC star formation rate is instructive. When all 
LMC data are included, we estimate the LMC star for¬ 


mation rate to be 0.3 M^:^ yr ^, assuming pjyrw = 1 (Ro- 
bitaille & Whitney 2Q1Q|) . This is within the error on 


the rate measured by [Harris & Zaritsky (2009), 0.4 ±q;2 
Mq (see their Figure 11), and given the magnetic field 
estimate performed in the previous section, we conclude 
that leptons make a nontrivial contribution to the LMC 
gamma-ray emission. 

In the case of 30 Dor, our s tar formation estimate is 
^2 times the rate measured by Harris & Zaritsky (2009) 
(see their Figure 17). If we take the star formation rate 
to be the measured value of 0 . 11 ±o;o 2 i^^ie lep¬ 

tonic fluxes estimated by our spectral fit both decrease by 
50%. The total estimated flux of 2.5 x 10“^ ph cm“^ s“^ 
is maintained given the anticorrelation between the nor¬ 
malizations of the leptonic and hadronic spectra. Taking 
this into account, we find that cosmic-ray electrons con¬ 
tribute ^45% of the 30 Dor emission. This wi ll affect the 
cos mic-ra y diffusion length scales reported by [Murphy et[ 
al. (2012), who assumed the 30 Dor gamma-ray emission 


was entirely due to hadronic collisions. 

We test the effect of the hydrogen density estimates 
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on the predicted leptonic contribution to the gamma-ray 
flux. The neutral and ionized hydrogen densities affect 
the electron energy los s rates due to bremsstrahlung and 
ionization (Equations |B1 Q}|bT2 ). As an alternative to 
= (2.0,0.0) cm“^ used in the analysis de¬ 
scribed above, we use (nui^nua ) = (0.06,0.06) cm~^. 
These densities were adopted by [Chakraborty fc Fields 


(2013) for their model of the inverse Compton gamma- 


ray spectrum in the Milky Way. Those authors consid¬ 
ered that electrons would not spend most of their time 
in the highest densities in the mid-plane of the Galactic 
disk but that magnetic fields cause electrons to migrate 
to regions of lower hydrogen density. 

When taking rim = = 0.06 cm“^, we find that 

the bremsstrahlung gamma-ray spectral model peaks at 
lower energies (^200 MeV). The best fit to the LMC 
gamma-ray spectrum using the new bremsstrahlung 
spectral model predicts a formation rate of 0.07 Mq yr“^, 
which is inconsistent with the rate of 0.4±o;2 mea¬ 


sured by Harris & Zaritsky (2009). However, if we fix 
the star formation rate parameter to 0.3 Mq yr“^ as de¬ 
termined by our previous fit , we find that the leptonic 
contribution to the LMC gamma-ray flux is still 30%. 

Our spatial analysis also suggests a strong leptonic 
component to the LMC gamma-ray flux. We found from 
our gas + radiation models that inverse Compton can ex¬ 
plain up to ^50% of the LMC gamma rays. Admittedly, 
this is an overestimate given that electron calorimetry, 
which was not considered for the spatial maps, would re¬ 
quire a corresponding bremsstrahlung component >50%. 
This imperfection aside, both our analysis methods agree 
that the the leptons contribute nontrivially to the LMC 
signal detected by Fermi. 

8.3. Diffusion of Cosmic-rays in the LMC 

We have convolved the star formation rate ma ps shown 
Figure!^ w ith the smoothing kernels used by |Mnrphy| 


et al.j (2012). We have the assembled a spatial model 
that includes the three gamma-ray emitting cosmic ray 
collision processes. This map takes the form: 

\J ^Ha(<^l'06.3 + CL2'f>12.b) * /^(^CRp) + 

a/ 7Ha(<^l'06.3 + O^2'012.5) * l^{lcRe) + 

^ ^160Mm(^lV^6.3 + ^2V^12.5) * /^(^CRe), 

(19) 

where the fluxes from each component are taken from 
the first row of Table [4l the and bi values are those 
given in Tables and |2fTor their re spective target model , 
and * is the convolution operator. Murphy et al. (2012) 
use two different smoothing kernels, n: the first is a two- 
dimensional gaussian, k{ 1) = e~^ , which models dif¬ 

fusion; the second is an exponential profile, k(V) = e~^^f 
which represents diffusion modified by cosmic-ray energy 
loss and escape. 


In their study, Murphy et al. (2012) found the best 
fit scale lengths for cosmic-ray electrons to be 0.1 and 
0.2 kpc for the exponential and gaussian kernels respec¬ 
tively. The scale lengths for protons were found to be 
0.2 and 0.45 kpc. Keepi ng the ratio /cRn/^CR e fixed to 
the value determined by Murphy et al. ([2012 ), we vary 
the length scales by factors of two, apply them to the 


model described by Equation (19), and run each model 
through the binned likelihood analysis discussed above. 
Figure shows the rankings of various models using 
diffusion smoothing. 



Fig. ff.— Rank of diffusion models as determined by the like¬ 
lihood statistic. Triangles indicate rankings based on the original 
dataset. Error bars indicate ±fcr spread on mean rankings deter¬ 
mined using ten boo tstrapped sarnples. l^/^ indicates the best fit 
length scale found by|Murphy et al.|j20f2|), and /^exp and ^Cgau refer 
to the exponential arid tiaussian smootJiing kernels respectively. 


All models that incorporate smoothing kernels outrank 
their unsmoothed counterparts. We find that do ubling 
and quadrupling the scaling lengths reported by |Mnr-| 
phy et al. (|2012) yield the highest ranked models. The 
exponential kernel with (/cRp, ^CRe) = (0.8,0.4) kpc 
is consistently ranked highest among both the original 
data set and the bootstrappe d samples. Thi s is at least 
in part due to the fact that Murphy et al. (2012) find 
^CRe ~ 0.5/cRp- Those authors determined only /cRp 
from gamma-ray data, but we are using both /cRp and 
^CRe to model the emission measured by Fermi. 

We note that this is a rather crude estimate of diffusion 
across the entire LMC. Spatial variation of the diffusion 
scale lengths can be caused by the magnetic structure of 
the LMC and may play an important role in our results; 
however, our current analysis cannot address this issue. 
The complexity of cosmic-ray propagation on the scale of 
an entire galaxy would be best studied using kinematic 
simulations such as GALPROP or idea l ly hydrodynam- 


ical s imulations (e.g. Booth et al. 2013 Salem & Bryan 
2014), which is beyond the scope of this article. 


8.4. 30 Dor as a Point Source and Additional LMC 
Model Parameters 

The previous analyses of the LMC Fermi data have 
been careful to consider the possibility that the gamma- 
ray emission from 30 Dor may not originate from cosmic- 
ray interactions, but rather from objects known to pro¬ 
duce point-like emission in the gamma-ray sky. Specifi¬ 
cally, the presence of pulsars and an X-ray background 
source are the causes for concern. Recently, HESS 
has observed one of these pulsars, PSR J0537-6910 and 
its associated wind nebula N 157B, at TeV energies 
(|H. E. S. S. Collaboration et al.||2015|). The spectrum 
they fit to the M 157B emission (their Figure 4) suggest 
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a flux of about 2 x 10“^^ erg cm“^ at 20 GeV, which 
is comparable to our 30 Dor flux estimate at the same 
wavelength (see Figure 10). However, at 4 GeV, their 
spectrum suggests the nux from N 157B should be an 
order of magnitude below the 30 Dor flux measured by 
Fermi. This pulsar cannot account for the bulk of the 30 
Dor gamma-ray emission at Fermi wavelengths. 

In a talk presented at the Fifth International Fermi 
Symposiuujj Pierrick Martin presented a spatial model 
used by the Fermi group for their recent analysis of 
gamma-ray emission from the LMC (publication in 
preparation, we thank our referee for pointing us to these 
slides). Their model includes two point sources that he 
outside the 30 Dor region. The first is supernova rem¬ 
nant N 132D, the type of object that has motivated our 
spatial maps because we expect it to be accelerating cos¬ 
mic rays. The location of N 132D already corresponds 
to a region of heightened star formation during the 12.5 
Myr epoch (Figure®. The second source also coincides 
with a region of enhanced star formation north of 30 
Dor, but has not been positively identified as a super¬ 
nova remnant. Based on the location indicated on slide 
5 of the talk referenced above, we have tested the effect 
of modeling this additional point source at the position 
of 2XMM J053444.6-673856. We find the rankings of 
our gas models remain unchanged, but the LMC flux es¬ 
timates decrease by ^8%. For readers familiar with the 
Fermi test statistic, the TS for this source is ^90, i.e., the 
log-likelihood increases by ^45 when we add this point 
source to the model. 

We have also tested the effect of adding an additional 
spectral parameter to our gas templates. The extra pa¬ 
rameter is defined such that the spectrum takes the shape 
of a parabola in log-space, and Equation ^ takes the 
form 


Ij{x,y,E) = y^aiMi{x,y)E 


-b^-Ci ln(_E) 


( 20 ) 


As was the case for the additional point source, we find 
the rankings of the gas models remain unchanged by the 
addition of the q parameter. The log-likelihoods increase 
by r\j 20 (TS^ 40) and the flux estimates are unchanged 
within error. 


9. CONCLUSIONS 

The LMC is the brightest and most spatially extended 
source of diffuse gamma-ray emission outside the Milky 
Way. It provides a unique laboratory for the study of 
cosmic-ray physics on global scales. The Fermi LMC ob¬ 
servations are thus complementary with the Milky Way 
diffuse emission, where smaller structures and even in¬ 
dividual supernova remnants can be resolved but where 
distances to sites of cosmic-ray acceleration may not be 
well constrained. 

We have performed both spatial and spectral analy¬ 
ses of the 5.5 years gamma-ray data from the LMC col¬ 
lected by the Fermi Gamma-ray Spaee Teleseope. Our 
spatial analysis proceeded by modeling the gamma-ray 
data with observed maps of the LMC radiation field and 
ISM components that are expected to contribute to the 

® http://fermi.gsfc.nasa.gov/science/mtgs/symposia/2014/ 
program/05_Martin.pdf 


gamma-ray emission, namely, star formation and non- 
thermal radio emission (cosmic-ray flux), ionized, neu¬ 
tral, and molecular hydrogen (target protons for tt^ de¬ 
cay and bremsstrahlung channels), and dust and infrared 
stellar emission (target photons for inverse Compton 
channel). We quantitatively distinguish between these 
models by ranking them by their likelihood statistic com¬ 
puted by the binned maximum likelihood estimator im¬ 
plemented in the Fermi Science Tools. We have tested 
the robustness of these rankings by applying bootstrap 
resampling to the LMC gamma-ray data set. 

From our spatial analysis, we conclude that cosmic rays 
remain relatively concentrated near their sites of accel¬ 
eration before losing energy via gamma-ray production. 
We find the gamma-ray spectral index to range between 
2.2 and 2.4 across all of our spatial models. In the case 
where we simultaneously fit gas and radiation models, we 
find inverse Compton contributes significantly to gamma- 
ray flux. A model for inverse Compton emission had not 
been considered by previous studies of the Fermi LMC 
data. 

We perform our sp ectra l analysis by first choosing a 
spatial model (see § 5.4 for model details) solely for 


the purpose of estimating LMC gamma-ray fluxes in 
six logarithmically spaced energy bins ranging from 200 
MeV to 20 GeV. To this gamma-ray spe ctrum, we fit 


a spectral m odel which includes tt^ decay ( Pfrommer fc 


EnfilinPOOd), inverse Compton scattering (Chakraborty 


.. 3), an d brems strahlung (formulated in Ap¬ 
pendix [b| (see Eqs. 14, and [D- Our spectral model 


fit uses least squares, which incorporates errors in the 
fluxes, to estimate the model parameters. 

Erom our spectral analysis, we find that the cosmic-ray 
proton power-law spectral index is ^2.4 to within 13% 
error. This value is consistent across all three sets of 
gamma-ray data, full LMC, 30 Dor removed, and 30 Dor 
isolated. The consistency leads us to conclude that vari¬ 
ation in the shape of the cosmic-ray spectrum is not the 
cause of the spatial variation in gamma-ray flux across 
the the LMC. Our spectral fits suggests that leptonic 
processes contribute a significant fraction of the LMC 
gamma-ray flux, which has not been considered by pre¬ 
vious studies. Assuming cosmic-ray energy equipartition 
with magnetic fields, we compute an LMC magnetic field 
strength of 2.8 /iG, wh ich is in good agreeme nt with the 
measurement made by Gaensler et al. (2005). _ 

We have appli ed the smoothing kernels used by |Mur 


|phy_et al. (2012) to the star formation rate maps in Eig- 
ure|^to rnodel cosmic-ray diffusion in the LMC. We find 
that d oubling and quadruplin g the scale lengths reported 
by the Murphy et al. (2012) result in the highest rank¬ 
ing models. We expect this is partly due to the fact 
that those authors determined only the cosmic-ray pro¬ 
ton scale length from the Fermi data, whereas we have 
fit the data using models for both electron and proton 
diffusion. 

We thank the anonymous referee for suggestions that 
vastly improved our paper. We thank Sui Ann Mao and 
Robert Brunner for a useful discussions and Keith Bech- 
tol for helpful recommendations about performing data 
reduction and analysis using the Fermi Science Tools 
software. This work is supported by NASA under the 
Fermi Guest Investigator Program (NNXIIAOI8G and 
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APPENDIX 

A. NON-THERMAL RADIO MAP 

To disting uish the thermal and nonT hermal components of the 1.4 GHz radio emission in the LMG, we use the method 
presented by Tabatabaei et al. (2007). In essence, this method constructs a model of the thermal radio emission from 
an extinction-corrected image of a galaxy’s Ha emission. The extinction correction is derived by estimating the dust 
optica l depth from the 160/im flux density and the dust temperature. The method is described in detail in Sections 
3-6 of Tabatabaei et al. (2007). Here we summarize the key equations, and discuss the values of specific parameters 
that we have adopted for the decomposition of the 1.4 GHz LMG m osaic. 


To estimate the dust temperature, we use the map constructed by Bernard et al. (2008). From this map, we obtain 
a map of the dust optical depth at 160 fim (rieo) according to: 


7i60 = Bieo{T) [1 - exp(ri6o)], 


(Al) 


where Iiqq is the flux density at 160 /im, and Biqq{T) is the value of the Planck function at 160 /im for dust temperature 
T. The dust optical depth at 160 /im is typically small, reaching a maximum value of rieo ^ 2.5 x 10“^ for several 
locations within the chain of molecular clouds south of 30 Doradus. 


Tabatabaei et al. (2007) use the extinction curve of the standard model for dust in the diffuse ISM to convert between 
ri6o and the dust optical depth at the wavelength of Ha emission adopting th^ ^ 2200ri6o- Since the large-grain 
populations in the LMG and Milky Way are likely to have similar optical properties (compare, for example, Figures 
7 and 8 of [Pgi|[T^ , we use the same conversion factor here. If the sources of Ha emission were located behind the 
galaxy, then the observed Ha emission /hq would be related to the intrinsic (i.e. extinction-free) emission /Ha,o via 


= ^Ha,0 exp(-rHa) 


(A2) 


We use the SHASSA map described in Section to estimate lua at each map pixel. In general, Ha sources lie within 
the galaxy, so r^a provides only an upper limit to the attenuation. The effective optical depth is Teff = /d n 
where /d G [0,1] is the dust-screening factor that represents the relative geometry of the Ha emission and the dust 
that contributes to th e extinction. If the dust, Ha sources, and DIG are well mixed, /d = 0.5. For the Milky Way, 
Dickinson et al. ([2003) find /d = 0.33 ± 0.15, indicating that Ha emission has a smaller vertical scale height than the 
dust. 4'he porous appearance of the HI emission in the LMG suggests that the ISM transparency in the LMG may 
be greater than in the Milky Way. We adopted /d = 0.1, which corresponds to a mean extinction in the LMG of 
Ana = 0.2 mag. 

Having obta ined an estimate for the intrinsic Ha flux density, lua = ^Ha,o exp(—Teff), we use equation 9 of Vails 
(1998) to estimate the emission measure EM: 


Gabaud 


^Ha,0 = 9.41 X 10 


Brp-1.017 


10 


0.029 

■ ^e4 EM. 


(A3) 


In this equation, Te 4 is the electron temperature Tg in units of 10^ K, and the expression is determined assuming Gase 
B recombination (i.e., each Lyman lin e photon is resonantly scattered many times). For twelve HII regions in the 


LMG, Vermeij & van der Hulst (2002) derived a mean electron temperature of 10,000 K. Individual measurements 
varied between 8000 and 16,000 K, depending on the location of the HII region and emission line that they used in 
their analysis. The electron temperature of the D IG in the LMG is not well determ ined; typical estimates in the 
Milky Way range between 8000 and 10,000K (e.g., Reynolds 1985 Alves et al. 2010). For our LMG decomposition. 


we adopted Tg = 8000 K. 

The op tical depth of the radio continuum emission Tg at frequency u is related to the emission measure derived from 
Equation |A3| by 

Tc = 8.235 X lO-^areiVGHz (1 + 0.08)£;M, (A4) 

where z^ghz is the observed freq uency expressed in GH z, and a is a correction factor that is approximately equal to 
unity at 1.4 GHz (see Table 3 of [Dickinson et al.||2003|). The predicted brightness temperature of the free-free radio 
continuum emission Tb is then simply 

Tb=Tg(l-exp(rg)). (A5) 

To obtain a map of the non-thermal radio continuum emission at 1.4 GHz, we subtract this model of the thermal radio 
emission from the median-filtered 1.4 GHz continuum map. The resulting maps of the LMG’s non-thermal 1.4 GHz 
radio continuum emission are presented in Figure The integrated non-thermal flux density is 265 Jy, corresponding 
to a global thermal fraction of 27%. 

The non-thermal map in Figure exhibits a mixture of diffuse emission and high-brightness features. On one 
hand, this finding could indicate that our thermal/non-thermal decomposition is failing for HII regions, i.e. that we 
are underestimating the Ha extinction and incorrectly identifying some of the bright thermal 1.4 GHz emission as 
being of non-thermal origin. On the other hand, some spatial coupling between the brightest sources of thermal and 
non-thermal 1.4 GHz emission would be expected since the massive stars that ionize the hydrogen gas in HII regions 
should evolve quickly to become supernova remnants. Overall, the brightest sources of non-t hermal radio emission 
demonstrate a good correspondence to the positions of confirmed LMG supernova remnants (Badenes et al. 2010), 
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Fig. 12.— Non-thermal 1.4 GHz map of the LMC used for our diffuse cosmic-ray electron models. See|Hughes et al.|(|2007|) for details 
about the original radio data. 


which should be strong synchrotron emitters. We find that < 20% of the LMC’s non-thermal emission arises from 
regions where the non-thermal radio flux density is greater than 0.2 mJy beam“^. A similar fraction is obtained if we 
calculate the non-thermal flux density above an Ha surface brightness threshold of 100 Rayleigh. The majority of the 
non-thermal radio emission in the LMC therefore arises from a diffuse, low surface brightness component that is not 
directly connected to HII regions. 

The method that we have used to separate the thermal and non-thermal components of the radio continuum emission 
has a number of limitations. In particular, we assume that the values of the dust-screening factor and the electron 
temperature across the LMC are constant even though both quantities should vary with interstellar environment. As 
noted above, Tg may be systematically lower in the HII regions tha n in the DIG. There is also empirical evidence for Tg 
fluctuations within individual HII regions (e.g. Tsamis et al.||2QQ4 ). The dust screening factor is probably not uniform 
either: the dust distribution may be more clurnpy in star-forming regions than in the diffuse ISM for example. For a 
constant dust mass, a clumpy distribution is more porous than a uniform layer, so the effective Ha attenuation will be 
higher in regions where the dust and gas are better mixed. A more sophisticated model would allow /d and Tg to vary 
across the LMC, but such a model would be difficult to justify without further independent constraints on /d and Tg. 


B. BREMSSTRAHLUNG SPECTRUM 

To derive Equation 0, we begin with the gamma-ray emissivity, which can be shown to equal 


_ . dN^ Aa ^ /-.ooN f 

QbremiE^) = ln(183)nH 


dE‘ 


Ae{E') 




A^a 

= —(Jt ln(I83)nH 


-^ 7 ) 


(Bl) 


Stecker 


from Equatio ns (8) to (36) of Stecker R97l|). We have taken the ultrarelativistic limit a ^rUgC^ (Equation 8 -28 

th' 


1971), and we have assumed the interstellar material is dominated by hydrogen, i.e., Z ^ 1. In Equation |B I 


0g is the cosmic-ray electron flux spectrum, a is the fine structure constant, and ctt is the Thomson scattering cross 
section. 

We adopt the cosmic-ray electron spectrum and the one-zone model used by Chakraborty & Eields (2013) so that 
we may consistently combine our bremsstrahlung gamma-ray spectrum with their inverse Compton spectrum. The 
cosmic-ray electron injection spectrum is given by 


qe{Ee) 


dN, 


dV dt dEe 


= TsNA4(Eg)E 


-1 


(B2) 


where Tsn is the supernova rate, which is proportional to the star formation rate -0, A4 is the number of electrons per 
unit ener gy accelerated by ea ch supernova event, and V is the galactic volume. A7g is taken to be a broken power law 
following Strong et al. (2010) 

Ee< Eh =4 GeV 

P-2.25 IP ^ p, _ /I 


Afe{Ee) OC 


E: 


where E^. is the cosmic-ray electron energy and Ei, is the power law break energy. We also assume a hard cutoff 
at Eg = 2 TeV. Taking the LMC as an electron calorimeter and the injection spectrum to be in steady state, the 
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cosmic-ray electron flux is given by 

4>e{Ee) = 


dN,. 


dA dt dEe 


[ dE'MK) = c 

JE. 


Ee) 


(B4) 


h{E,)JE^ h{E,) 

where we have assumed Ef. ^ rtieC?, and b{E^) incorporates all the energy-loss mechanisms available for cosmic-ray 
electrons i e 

b{Eg) = bic{Eg) + bsynciEe) + ^brem(-E’e) + ^ion(-Be)- (B5) 

For the injection spectrum given by Eqs. ( |B2| ) and ( |B3| ), 

p-0.s^ 


qe{> Ee) OC 


(a 


25(£;-» s - El 


r-i.25\ 


Chakraborty fc Fields| ( 2Q13[ ) give the synchrotron energy-loss rate as 
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where B is the magnetic field strength. In the Thomson limit where electron energies are low enough that the Klein- 
Neshina correction to the inverse Compton cross section is unimportant, 


^ic(^e) = 2.5 


b^ISRF 


1.1 eV cm“^ 


4/iG 

B 


^^sync 
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(B8) 


where Uisrf is the energy density of the ISRF. The bremsstrahlung energy loss rate is given as the sum of two 
components, 

^brem(-F'e) = 5brem,i(F'e) T 5brem,n(F'e), (B9) 

where 
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is the contribution from cosmic-ray electron interactions with ionized hydrogen and 

6brem.n(^e) = 7.3 X lO’^" GcV s" 


/ nm \ 
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(Bll) 

i s the contribution from neutral hydrogen. We use the ionization energy loss term given by [Ginzburg & Syrovatskii| 
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B.l. Milky Way Bremsstrahlung 


As did Chakraborty & Fields (2013) f or the inverse Compt on spectrum, we normalize the bremsstrahlung spectrum 
to that of t he plain dihnsion m odel of [Strong et al.| (|2010|) (see solid cyan curve of right panel of Figure 1). The 
spectrum of Strong et al. (2010) is given in terms of specihc luminosity 


Fbr 


liEj) — J Qhrem.{Ee)dV — ^brem(-F'e) J dV^ 


(B13) 


where the seco nd equality incorporates the one-zone model a ssum ption. Notice here that the volume integral in 
Equation (B13) cancels with the factor of V~ ^ in E q uatio n (B2 ), wh ich means the bremsstrahlung luminosity is 
independent of galactic volume. Erom Eqs. (|B4|), (B6), (Bl), arm^BlS) 
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(B14) 

where all energies are expressed in GeV and A* is the constant to be determined by normalizing to the [Strong et al.| 
(2010) bremsstrahlung spectrum and is proportional to the electron-injection rate, which is proportional to the star 
fo rmation rate. Again, Eh = 4 G eV and Emax = 2 TeV. Eor the Milky Way energy loss rates, we adopt the same values 


as 

the 

0-1 


Chakraborty fc Eields (2013): Visrf = lA eV cm ^,5 = 4 /iG, and n m = nua = O.Snp = 0.06 cm Based on 
Strong et al. (2010|) bremsstrahlung spectrum, we normalize Equation (B14) such that E^I/brem(^e) = 10^^ GeV^ 
GeV“^ at e\ = 200 MeV, which results in A* = 5.93 x 10^^ s“^. 
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B.2. LMC Bremsstrahlung 

The relative importances of the electron energy-loss mechanisms are different in the LMC from those in the Milky 
Way. In particula r, bremsstrahlung losses are more important given the higher neutral-hydrogen density nni = 2 cm“^ 
(Kim et ah 2003). Inverse Compton losses are less important given t he lower ISRF energy density I/isrf = 0.57 eV 
cm“'^, calculated from the LMC SED presented by Israel et al. (2010) where we added the CMB contribution us ing a 
blackbody spectrum with T = 2.73 K. For the magnetic held strength, we’ve used B = A /iC (Gaensler et al. ]200^ . We 
also assume the ionized-hydrogen density is negligible compared with the HI density. We use these quantities in the 
energy-loss-rate equations, and we parameterize the normalization of LMC bremsstrahlung spectrum as V^lmc/V^mw 
resulting in Equation ( pT| ). 
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